%*********************************Solving Generalized Moving Peaks Benchmark (GMPB)_SIMPLIFIED VERSION (only one sub-function, i.e., non-modular) using mQSO******
%Author: Danial Yazdani
%Last Edited: June 03, 2021
%
% ------------
% Reference:
% ------------
%
%  D. Yazdani et al.,
%            "Benchmarking Continuous Dynamic Optimization: Survey and Generalized Test Suite,"
%            IEEE Transactions on Cybernetics (2020).

%  T. Blackwell and J. Branke,
%            "Multiswarms, exclusion, and anti-convergence in dynamic environments"
%            IEEE Transactions on Evolutionary Computation (2006).
% ------------
% Notification:
% ------------
% It is assumed that the environmental changes are VISIBLE, % therefore,
% mQSO is informed about changes (i.e., mQSO does not need to detect
% environmental changes). Also note that mQSO does not access to a prior knowledge
% about the shift severity value. The shift severity is learned in this code.
%
% --------
% License:
% --------
% This program is to be used under the terms of the GNU General Public License
% (http://www.gnu.org/copyleft/gpl.html).
% Author: Danial Yazdani
% e-mail: danial.yazdani AT gmail dot com
%         danial.yazdani AT yahoo dot com
% Copyright notice: (c) 2021 Danial Yazdani
%**************************************************************************************************
clear all;close all;clc;
%% Input variables
%********Figures and visualizations*************
OfflineErrorFig                          = 0;   % set to 1 to draw offline error over time plot, 0 otherwise
CurrentErrorFig                          = 1;   % set to 1 to draw Current error plot, 0 otherwise
VisualizationErrorOverOptimization       = 0;   %set to 1 to draw Current error and offline error plot in real-time, 0 otherwise
VisualizationIndividualsOverOptimization = 0;   %set to 1 to see the contour of the landscape and individuals in real-time, 0 otherwise. Note that this only works in 2-dimensional space
%********Benchmark parameters and Run number***
PeakNumber                               = 10;  %The default value is 10
ChangeFrequency                          = 5000;%The default value is 5000
Dimension                                = 5;   %The default value is 10
ShiftSeverity                            = 1;   %The default value is 1
EnvironmentNumber                        = 100;  %The default value is 100
RunNumber                                = 5;   %It should be set to 31
%%
BestErrorBeforeChange = NaN(1,RunNumber);
OfflineError = NaN(1,RunNumber);
OfflinePerformance = NaN(1,RunNumber);
CurrentError = NaN (RunNumber,ChangeFrequency*EnvironmentNumber);
for RunCounter=1 : RunNumber
    rng(RunCounter);%This random seed setting is used to initialize the Problem-This must be identical for all peer algorithms to have a fair comparison.
    Problem = BenchmarkGenerator(PeakNumber,ChangeFrequency,Dimension,ShiftSeverity,EnvironmentNumber);
    rng('shuffle');%Set a random seed for the optimizer based on the system clock
    %% Initialiing Optimizer
    clear Optimizer;
    Results = NaN(2,Problem.EnvironmentNumber);
    Optimizer.Dimension = Problem.Dimension;
    Optimizer.PopulationSize = 5;
    Optimizer.MaxCoordinate   = Problem.MaxCoordinate;
    Optimizer.MinCoordinate = Problem.MinCoordinate;
    Optimizer.DiversityPlus = 1;
    Optimizer.x = 0.729843788;
    Optimizer.c1 = 2.05;
    Optimizer.c2 = 2.05;
    Optimizer.ShiftSeverity = 1;%initial shift severity
    Optimizer.QuantumRadius = Optimizer.ShiftSeverity;
    Optimizer.QuantumNumber = 5;
    Optimizer.SwarmNumber = 10;
    Optimizer.ExclusionLimit = 0.5 * ((Optimizer.MaxCoordinate-Optimizer.MinCoordinate) / ((Optimizer.SwarmNumber) ^ (1 / Optimizer.Dimension)));
    Optimizer.ConvergenceLimit = Optimizer.ExclusionLimit;
    for ii=1 : Optimizer.SwarmNumber
        [Optimizer.pop(ii),Problem] = InitializingOptimizer(Optimizer.Dimension,Optimizer.MinCoordinate,Optimizer.MaxCoordinate,Optimizer.PopulationSize,Problem);
    end
    VisualizationFlag=0;
    %% main loop
    while 1
        %% Visualization
        if (VisualizationIndividualsOverOptimization==1 && Dimension == 2) && VisualizationErrorOverOptimization == 1
            ax1 = subplot(1,2,1);
            ax2 = subplot(1,2,2);
        elseif (VisualizationIndividualsOverOptimization==1 && Dimension == 2)
            ax1 = subplot(1,1,1);
        elseif VisualizationErrorOverOptimization == 1
            ax2 = subplot(1,1,1);
        end
        if VisualizationIndividualsOverOptimization==1 && Dimension == 2
            if VisualizationFlag==0
                VisualizationFlag=1;
                T = Problem.MinCoordinate : ( Problem.MaxCoordinate-Problem.MinCoordinate)/100 :  Problem.MaxCoordinate;
                L=length(T);
                F=zeros(L);
                for i=1:L
                    for j=1:L
                        F(i,j) = DUMMYfitness([T(i), T(j)],Problem);
                    end
                end
            end
            contour(ax1,T,T,F,25);
            colormap cool;
            xlabel(ax1,'x_1')
            ylabel(ax1,'x_2')
            grid(ax1,'on');
            hold(ax1,'on');
            for ii=1 : PeakNumber
                if Problem.PeakVisibility(Problem.Environmentcounter,ii)==1
                    if ii == Problem.OptimumID(Problem.Environmentcounter)
                        plot(ax1,Problem.PeaksPosition(ii,2,Problem.Environmentcounter),Problem.PeaksPosition(ii,1,Problem.Environmentcounter),'p','markersize',15,'markerfacecolor','none','MarkerEdgeColor','k','LineWidth',1.5);
                    else
                        plot(ax1,Problem.PeaksPosition(ii,2,Problem.Environmentcounter),Problem.PeaksPosition(ii,1,Problem.Environmentcounter),'o','markersize',15,'markerfacecolor','none','MarkerEdgeColor','k','LineWidth',1.5);
                    end
                end
            end
            for ii=1 : Optimizer.SwarmNumber
                for jj=1 : Optimizer.PopulationSize
                    plot(ax1,Optimizer.pop(ii).X(jj,2),Optimizer.pop(ii).X(jj,1),'o','markersize',7,'markerfacecolor','g','MarkerEdgeColor','none');
                    plot(ax1,Optimizer.pop(ii).BestPosition(2),Optimizer.pop(ii).BestPosition(1),'p','markersize',10,'markerfacecolor','g','LineWidth',1.2);
                end
            end
            hold(ax1,'off');
        end
        if VisualizationErrorOverOptimization == 1
            semilogy(ax2,Problem.CurrentError,'r','DisplayName','Current Error');
            xlabel(ax2,'Fitness Evaluation');
            ylabel(ax2,'Current Error');
            xlim(ax2,[0 Problem.MaxEvals])
            ylim(ax2,[0 1000])
            legend;
            grid(ax2,'on');
        end
        if (VisualizationIndividualsOverOptimization==1 && Dimension == 2) && VisualizationErrorOverOptimization == 1
            set(gcf, 'Position',  [100, 100, 1400, 600])
            drawnow;
        elseif (VisualizationIndividualsOverOptimization==1 && Dimension == 2)  || VisualizationErrorOverOptimization == 1
            drawnow;
        end
        %% Optimization
        [Optimizer,Problem] = Optimization(Optimizer,Problem);
        if Problem.RecentChange == 1%When an environmental change has happened
            Problem.RecentChange = 0;
            [Optimizer,Problem] = Reaction(Optimizer,Problem);
            VisualizationFlag = 0;
            clc; disp(['Run number: ',num2str(RunCounter),'   Environment number: ',num2str(Problem.Environmentcounter)]);
        end
        if  Problem.FE >= Problem.MaxEvals%When termination criteria has been met
            break;
        end
    end
    BestErrorBeforeChange(1,RunCounter) = mean(Problem.Ebbc);
    OfflineError(1,RunCounter) = mean(Problem.CurrentError);
    OfflinePerformance(1,RunCounter) = mean(Problem.CurrentPerformance);
    CurrentError(RunCounter,:) = Problem.CurrentError;
end
E_bbc = [mean(BestErrorBeforeChange),median(BestErrorBeforeChange),std(BestErrorBeforeChange)/sqrt(RunNumber)];
E_o = [mean(OfflineError),median(OfflineError),std(OfflineError)/sqrt(RunNumber)];
P_o = [mean(OfflinePerformance),median(OfflinePerformance),std(OfflinePerformance)/sqrt(RunNumber)];
close;clc;
disp(['Offline error ==> ', ' Mean = ', num2str(E_o(1)), ', Median = ', num2str(E_o(2)), ', Standard Error = ', num2str(E_o(3))]);
disp(['Offline Performance ==> ', ' Mean = ', num2str(P_o(1)), ', Median = ', num2str(P_o(2)), ', Standard Error = ', num2str(P_o(3))]);
disp(['Average error before change ==> ', ' Mean = ', num2str(E_bbc(1)), ', Median = ', num2str(E_bbc(2)), ', Standard Error = ', num2str(E_bbc(3))]);
if CurrentErrorFig==1 || OfflineErrorFig==1
    FigureCurrentError = mean(CurrentError);
end
if CurrentErrorFig==1
    semilogy(FigureCurrentError,'r','DisplayName','Current Error');
    xlabel('Fitness Evaluation');
    ylabel('Error');
    set(gcf, 'Position',  [100, 100, 900, 400])
    legend;
    grid on;
end
if OfflineErrorFig==1
    if CurrentErrorFig==1
        hold on;
    end
    disp('Generating figure: please wait for calculating the offline error over time...');
    FigureOfflineError = NaN(size(FigureCurrentError));
    parfor ii=1 : length(FigureCurrentError)
        FigureOfflineError(ii)= mean(FigureCurrentError(1:ii));
    end
    semilogy(FigureOfflineError,'b','LineWidth',2,'DisplayName','Offline Error');
    xlabel('Fitness Evaluation');
    ylabel('Error');
    set(gcf, 'Position',  [100, 100, 900, 400])
    legend;
    grid on;
end